Consignes

Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :

La bioinfo c'est : <code>MERVEILLEUX</code>.

N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.

Ce rapport devra être mis à disposition et partagé avec nous sous sa forme compilée (html sous forme de Github Pages ou à défautPDF) dans votre dépôt github.

Les parties 1 à 4 seront notées pour l’évaluation du module 4, les parties 5 et 6 pour le module 5.

Introduction

Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :

Analyses

Organisation de votre espace de travail

mkdir ~/EVALUATION
cd ~/EVALUATION

Téléchargement des données brutes

  • Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]
module load sra-tools
fasterq-dump --version
# "fasterq-dump" version 2.10.3
srun --cpus-per-task 8 fasterq-dump -S -p SRR10390685 --outdir . --threads 8
gzip SRR10390685_1.fastq
gzip SRR10390685_2.fastq
  • Combien de reads sont présents dans les fichiers R1 et R2 ?
zcat SRR10390685_1.fastq.gz | echo $((`wc -l`/4))
zcat SRR10390685_2.fastq.gz | echo $((`wc -l`/4))
ou
module load seqkit
seqkit stat SRR10390685_*.fastq.gz
# file                    format  type   num_seqs        sum_len  min_len  avg_len  max_len
# SRR10390685_1.fastq.gz  FASTQ   DNA   7,066,055  1,056,334,498       35    149.5      151
# SRR10390685_2.fastq.gz  FASTQ   DNA   7,066,055  1,062,807,718      130    150.4      151

Les fichiers FASTQ contiennent 7 066 055 reads.

  • Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz
gzip -d GCF_000009045.1_ASM904v1_genomic.fna.gz
  • Quelle est la taille de ce génome ?
module load seqkit
seqkit version
# seqkit v0.14.0

zcat  GCF_000009045.1_ASM904v1_genomic.fna |seqkit stat 
# file                                  format  type  num_seqs    sum_len    min_len    avg_len    max_len
# GCF_000009045.1_ASM904v1_genomic.fna  FASTA   DNA          1  4,215,606  4,215,606  4,215,606  4,215,606

La taille de ce génome est de 4 215 606 paires de bases.

  • Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz
  • Combien de gènes sont connus pour ce génome ?
zgrep -v "^#" GCF_000009045.1_ASM904v1_genomic.gff.gz | awk '($3 == "gene")' |wc -l
# 4448

4 448 gènes sont recensés dans le fichier d’annotation.

Contrôle qualité

  • Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit
module load fastqc
fastqc --version
# FastQC v0.11.9

mkdir QC
srun --cpus-per-task 8 fastqc SRR10390685_1.fastq.gz -o QC/ -t 8
srun --cpus-per-task 8 fastqc SRR10390685_2.fastq.gz -o QC/ -t 8

module load multiqc
multiqc --version
# multiqc, version 1.9

srun multiqc -d QC -o QC
  • La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?

La qualité des bases me paraît … car … comme le montre …

Lien vers le rapport MulitQC

  • Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?

Oui/Non car …

  • Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?
(R1+R2) / Taille du génome

La profondeur de séquençage est de : 502.6898187 X.

Nettoyage des reads

Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp qui vous semblent adéquats et justifiez-les.

module load fastp
fastp --version
# fastp 0.20.0

mkdir FASTP
srun --cpus-per-task 8 fastp --in1 SRR10390685_1.fastq.gz --in2 SRR10390685_2.fastq.gz --out1 FASTP/SRR10390685_1.fastq.gz --out2 FASTP/SRR10390685_2.fastq.gz --html FASTP/fastp.html --thread 8 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json FASTP/fastp.json

seqkit stat FASTP/SRR10390685_[12].fastq.gz.

Les paramètres suivants ont été choisis :

Parametre Valeur Explication
–cut_mean_quality 30 pour un score moyen dans la fenêtre glissante > 30
–cut_window_size 8 pour une taille de fenêtre glissante de 8
–length_required 100 pour ne garder que les reads de taille > 100
–cut_tail pour faire partir la fenêtre de l’extrémité 3’ du read

Ces paramètres ont permis de conserver 6 777 048 reads pairés, soit une perte de 4.0900757 % des reads bruts.

Alignement des reads sur le génome de référence

Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [3] et samtools [4].

module load samtools
samtools --version
# samtools 1.10
# Using htslib 1.10.2

module load bwa
bwa
# Version: 0.7.17-r1188

srun bwa index GCF_000009045.1_ASM904v1_genomic.fna.gz
mkdir MAPPING
srun --cpus-per-task=4 bwa mem GCF_000009045.1_ASM904v1_genomic.fna FASTP/SRR10390685_1.fastq.gz FASTP/SRR10390685_2.fastq.gz -t 3 | samtools view -hbS - > MAPPING/SRR10390685.bam
srun samtools flagstat MAPPING/SRR10390685.bam > MAPPING/SRR10390685.bam.flagstat
srun samtools sort MAPPING/SRR10390685.bam -o MAPPING/SRR10390685_sorted.bam
srun samtools index MAPPING/SRR10390685_sorted.bam
  • Combien de reads ne sont pas mappés ?
samtools view -f 4 -c MAPPING/SRR10390685.bam
# 744540

744 540 reads ne sont pas mappés.

Croisement de données

Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [5]:

module load bedtools
bedtools --version
# bedtools v2.29.2

grep trmNF GCF_000009045.1_ASM904v1_genomic.gff | awk '$3=="gene"' > trmNF.gff3
bedtools intersect -a MAPPING/SRR10390685_sorted.bam -b trmNF.gff3.gz -f 0.5 > SRR10390685_on_trmNF.bam
samtools view -c SRR10390685_on_trmNF.bam

2 801 reads chevauchent le gène d’intérêt.

Visualisation

Utilisez IGV [6] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.

IGV : couverture du gène trmNF

References

1. toolkit NS. NCBI sra toolkit. NCBI, GitHub repository. 2019.

2. Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

3. Li H. Aligning sequence reads, clone sequences and assembly contigs with bwa-mem. arXiv preprint arXiv:13033997. 2013.

4. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and samtools. Bioinformatics. 2009;25:2078–9.

5. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.

6. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative genomics viewer (igv): High-performance genomics data visualization and exploration. Briefings in bioinformatics. 2013;14:178–92.

 

A work by Migale Bioinformatics Facility

https://migale.inrae.fr

Our two affiliations to cite us:

Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France

Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France